Demonstration: Compare Density Models¶
This notebook demonstrates how to compare the GEODYN output for different density models
First, load the data using the reading functions.
[1]:
%load_ext autoreload
%autoreload 2
import pandas as pd
import numpy as np
##################################
# INPUT PARAMETERS:
##################################
sat_file = 'starlette'
arc1 = '030914_2wk'
grav_id ='goco05s'
local_path = '/data/analysis/starlette_analysis/'
SAT_ID = 7501001
##################################
# PATH TO DENSITY MODEL RUN of Choice:
##################################
msis86_model = 'msis86'
path_to_msis86 = '/data/runs_geodyn/st/results/'+ grav_id+'_'+msis86_model+'/'
dtm87_model = 'dtm87'
path_to_dtm87 = '/data/runs_geodyn/st/results/'+ grav_id+'_'+dtm87_model+'/'
jaachia71_model = 'jaachia71'
path_to_jaachia71 = '/data/runs_geodyn/st/results/'+ grav_id+'_'+jaachia71_model+'/'
import sys
sys.path.insert(0, '/data/analysis/util_funcs/py_starlette/')
from a_ReadStarlette import ReadStarlette
(AdjustedParams_msis86,
Trajectory_msis86,
Density_msis86,
Resids_msis86) = ReadStarlette(arc1,
sat_file,
grav_id,
local_path,
path_to_msis86)
(AdjustedParams_dtm87,
Trajectory_dtm87,
Density_dtm87,
Resid_dtm87s) = ReadStarlette(arc1,
sat_file,
grav_id,
local_path,
path_to_dtm87)
(AdjustedParams_jaachia,
Trajectory_jaachia,
Density_jaachia,
Resids_jaachia) = ReadStarlette(arc1,
sat_file,
grav_id,
local_path,
path_to_jaachia71)
The base file name for this arc is: st030914_2wk.goco05s
File exists: iieout
/data/runs_geodyn/st/results/goco05s_msis86/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
/data/runs_geodyn/st/results/goco05s_msis86/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
/data/runs_geodyn/st/results/goco05s_msis86/DENSITY/st030914_2wk.goco05s
Loading data...
Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded
The base file name for this arc is: st030914_2wk.goco05s
File exists: iieout
/data/runs_geodyn/st/results/goco05s_dtm87/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
/data/runs_geodyn/st/results/goco05s_dtm87/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
/data/runs_geodyn/st/results/goco05s_dtm87/DENSITY/st030914_2wk.goco05s
Loading data...
Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded
The base file name for this arc is: st030914_2wk.goco05s
File exists: iieout
/data/runs_geodyn/st/results/goco05s_jaachia71/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
/data/runs_geodyn/st/results/goco05s_jaachia71/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
/data/runs_geodyn/st/results/goco05s_jaachia71/DENSITY/st030914_2wk.goco05s
Loading data...
Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded
Drag Coefficient Comparison¶
[7]:
labels = list(AdjustedParams_msis86[5]['0CD'].keys())
val_list = []
for i in AdjustedParams_msis86[5]['0CD'].keys():
val_list.append(AdjustedParams_msis86[5]['0CD'][i]['CURRENT_VALUE'])
[8]:
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import rc
import matplotlib as mpl
from matplotlib import patches as mpatches
mpl.rcParams['lines.markersize'] = 10
plt.rcParams['axes.grid'] = True
plt.rcParams['grid.alpha'] = 1
plt.rcParams['grid.color'] = "#cccccc"
plt.rcParams.update({'font.size': 16})
rc('font',**{'family':'sans-serif','sans-serif':['Calibri']})
plt.rcParams['axes.titlesize']='large'
plt.rcParams['axes.titlepad']= 3
rc('text', usetex=False)
plt.rcParams["legend.loc"] = 'best'
color86 = 'tab:blue'
color71 = 'tab:red'
color87 = 'tab:green'
fig, ( ax1) = plt.subplots(1, figsize=(7,6), sharex=False)
for i in AdjustedParams_msis86.keys():
ax1.set_title('Cd Output (Current Value) at each Iteration', y=1.1)
ax1.plot(i, AdjustedParams_msis86[i]['0CD'][labels[0]]['CURRENT_VALUE'],'.', color = color86, label = 'MSIS 86 T01')
ax1.plot(i, AdjustedParams_dtm87[i]['0CD'][labels[0]]['CURRENT_VALUE'],'x', color = color87, label = 'DTM T01')
ax1.plot(i, AdjustedParams_jaachia[i]['0CD'][labels[0]]['CURRENT_VALUE'],'v', color = color71, label = 'Jaachia T01')
ax1.set(ylabel= 'Cd [unitless]')
ax1.set(xlabel= 'Iteration Number')
msis_patch = mpatches.Patch(color=color86, label='MSIS')
dtm_patch = mpatches.Patch(color=color87, label='DTM')
jaachia_patch = mpatches.Patch(color=color71, label='Jaachia')
fig.axes[0].legend( handles=[msis_patch, dtm_patch, jaachia_patch],bbox_to_anchor=(0, 1.0, 1.0, .08), ncol=3, mode="expand", borderaxespad=0.)
fig.tight_layout( pad=1.0)
for ax in fig.axes:
plt.sca(ax)
plt.xticks(rotation=45)
findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans.
findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans.
The drag coefficient is time dependent¶
[9]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
from plotly.subplots import make_subplots
last_iter = 5
which_stat = 'CURRENT_VALUE' # 2 is the current val
labels = list(AdjustedParams_msis86[5]['0CD'].keys())
fig = make_subplots(
rows=1, cols=1,
subplot_titles=("Time Dependent Drag Coefficient",
))
val_list = []
for i in AdjustedParams_msis86[last_iter]['0CD'].keys():
val_list.append(AdjustedParams_msis86[last_iter]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
y=val_list,
name= 'MSIS 86',
mode='markers',
marker=dict(
size=10,),
), row=1, col=1,
)
val_list = []
for i in AdjustedParams_dtm87[last_iter]['0CD'].keys():
val_list.append(AdjustedParams_dtm87[last_iter]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
y=val_list,
name= 'DTM 87',
mode='markers',
marker=dict(
size=10,),
), row=1, col=1,
)
val_list = []
for i in AdjustedParams_jaachia[last_iter]['0CD'].keys():
val_list.append(AdjustedParams_jaachia[last_iter]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
y=val_list,
name= 'Jaachia 71',
mode='markers',
marker=dict(
size=10,),
), row=1, col=1,
)
fig.update_yaxes( title=r"$C_D (Unitless)$",exponentformat= 'power',row=1, col=1)
fig.update_xaxes( title="Iteration Number", row=1, col=1)
# fig.update_yaxes( yaxis_title=r"$\frac{m}{s^2}$", xaxis_title="Iteration Number",exponentformat= 'power',row=2, col=1)
# fig.update_layout(
# title="Along Track Accleration Amplitude vs Iteration", )
# fig.update_layout(
# yaxis_title=r"$\frac{m}{s^2}$",
# xaxis_title="Iteration Number" )
iplot(fig)
Compare density along the orbit¶
[10]:
data_nums = 1000
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
fig = go.Figure(data=[go.Scatter(x=Density_msis86['Date'][:data_nums],
y=Density_msis86['rho (kg/m**3)'][:data_nums],
name='MSIS 86',
mode='markers',
marker=dict(
size=2,
),
),
],
)
fig.add_trace(go.Scatter(x=Density_dtm87['Date'][:data_nums],
y=Density_dtm87['rho (kg/m**3)'][:data_nums],
name= 'DTM 87',
mode='markers',
marker=dict(
size=2,),
),
)
fig.add_trace(go.Scatter(x=Density_jaachia['Date'][:data_nums],
y=Density_jaachia['rho (kg/m**3)'][:data_nums],
name= 'Jaachia 71',
mode='markers',
marker=dict(
size=2,),
),
)
fig.update_layout(
title="Density along Starlette Orbit (for 4 minutes of first day)",
yaxis_title=r'$\frac{kg}{m^3}$',
xaxis_title="Date",
)
fig.update_yaxes(type="log", exponentformat= 'power',)
iplot(fig)
[11]:
data_nums = 15
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
fig = go.Figure(data=[go.Scatter(x=Density_msis86['Date'][::data_nums],
y=Density_msis86['rho (kg/m**3)'][::data_nums],
name = 'MSIS 86',
mode='markers',
marker=dict(
size=2,
),
),
],
)
fig.add_trace(go.Scatter(x=Density_dtm87['Date'][::data_nums],
y=Density_dtm87['rho (kg/m**3)'][::data_nums],
name= 'DTM 87',
mode='markers',
marker=dict(
size=2,),
),
)
fig.add_trace(go.Scatter(x=Density_jaachia['Date'][::data_nums],
y=Density_jaachia['rho (kg/m**3)'][::data_nums],
name= 'Jaachia 71',
mode='markers',
marker=dict(
size=2,),
),
)
fig.update_layout(
title="Density along Starlette Orbit (Every 15th datapoint)",
yaxis_title=r'$\frac{kg}{m^3}$',
xaxis_title="Date",
)
fig.update_yaxes(type="log", exponentformat= 'power',)
iplot(fig)
Compare residuals¶
[12]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
fig = go.Figure(data=[go.Scatter(x=pd.to_datetime(Resids_msis86['Date']),
y=Resids_msis86['Residual'].values.astype(float),
name ='MSIS 86' ,
mode='markers',
marker=dict(
size=3,
),
),
],
)
fig.add_trace(go.Scatter(x=pd.to_datetime(Resid_dtm87s['Date']),
y=Resid_dtm87s['Residual'].values.astype(float),
name= 'DTM 87',
mode='markers',
marker=dict(
size=3,),
),
)
fig.add_trace(go.Scatter(x=pd.to_datetime(Resids_jaachia['Date']),
y=Resids_jaachia['Residual'].values.astype(float),
name= 'Jaachia 71',
mode='markers',
marker=dict(
size=3,),
),
)
fig.update_layout(
title="Observation Residuals from Final Iteration",
yaxis_title="Residuals",
xaxis_title="Date",
)
iplot(fig)